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ABSTRACT 

The Tibet-Ill air shower array, consisting of 533 scintillation detectors, has been operating success- 
fully at Yangbajing in Tibet, China since 1999. Using the dataset collected by this array from 1999 
November through 2005 November, we obtained the energy spectrum of 7-rays from the Crab Nebula, 
expressed by a power law as (dJ/dE) = (2.09 ± 0.32) x 10~ 12 (£/3 TeV)- 2 - 96±014 cm- 2 s- 1 TeV _1 in 
the energy range of 1.7 to 40 TeV. This result is consistent with other independent 7-ray observations 
by imaging air Cherenkov telescopes. In this paper, we carefully checked and tuned the performance 
of the Tibet-Ill array using data on the moon's shadow in comparison with a detailed Monte Carlo 
simulation. The shadow is shifted to the west of the moon's apparent position as an effect of the 
geomagnetic field, although the extent of this displacement depends on the primary energy positively 
charged cosmic rays. This finding enables us to estimate the systematic error in determining the 
primary energy from its shower size. This error is estimated to be less than ±12% in our experiment. 
This energy scale estimation is the first attempt among cosmic-ray experiments at ground level. The 
systematic pointing error is also estimated to be smaller than 0?011. The deficit rate and position of 
the moon's shadow are shown to be very stable within a statistical error of ±6% year by year. This 
guarantees the long-term stability of point-like source observation with the Tibet-Ill array. These 
systematic errors are adequately taken into account in our study of the Crab Nebula. 
Subject headings: cosmic rays — gamma rays : observations — magnetic fields — Moon — pulsars : 
individual (Crab pulsar) — supernova remnants : individual (Crab Nebula) 
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1. INTRODUCTION 

The Crab Nebula is a standard source of radiation in 
the northern sky across a wide energy band, from radio 
to near 100 TeV 7-rays. It is well known that the multi- 
wavelength non-thermal energy spectrum is dominated 
by synchrotron radiation at energies lower than 1 GeV 
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and by the inverse-Compton scattering above 1 GeV 
(|De Jager et al. 199dlAtoyan fc Aharonian 199H . 

TeV 7-rays from the Crab Nebula were first 
clearly detected by the Whipple collaboration using 
an imaging air Ch erenkov telescope (1ACT) in 1989 
(jWeekes et al. 1989tl . Since then, IACT has become the 
standard telescope for high-energy 7-ray observations by 
virtue of its excellent angular resolution and efficiency. 
Many IACTs have been constructed and are operating 
around the world, detecting about 70 7-ray sources in 
total up to the present. On the other hand, air shower 
arrays have been constructed to search for 7-rays from 
point sources at high altitude. The merit of this tech- 
nique is that it can be operated for 24 hours every day, 
regardless of weather, with a wide field of view of about 
2 sr. The energy threshold of the 7-rays detected is 
higher than that of IACTs, say about 3 TeV at high 
altitude. Among these instruments, the Tibet AS7 Col- 
laboration achieved the first successful observation of the 
Crab Nebula at a multi-TeV region in 1999, using a so- 
called HD (high density) array in which 109 scintillation 
detectors were deployed at 7.5 m spacing lattice in tervals 
in an area of 5,175 m 2 (| Amenomori et al. 19 99a). The 
Milagro group also reported the detection of TeV 7-ray 
signa ls from the Crab N ebula using a water Cherenkov 
pool (| Atkins et al. 20031 ). 

Recently, IACTs have obtained updated informa- 
tion on the Crab Nebula. The HEGRA experi- 
ment has extended the nebula's energy spectrum up 
to 80 TeV with an approximate power-law shape, af- 
ter patient observation for almost 400 hours in total 
(| Aharonian et al. 200l) . In contrast, the MAGIC ex- 
periment, equipped with the world's largest tessellated 
reflector (17 m in diameter), has successfully observed 
the lower energy pa rt of the spectrum down to 77 GeV 
(| Albert et al. 2 008). The H.E.S.S. group examined the 
energy spectra of the Crab Nebula obtained from var- 
ious IACTs (Whipple, HEGRA, CAT and H.E.S.S.) 
to evaluate the system atic errors of these instruments 
(| Aharonian et al. 2 006). Their fluxes are well in agree- 
ment with one another within the statistical and system- 
atic errors at the moderate energy region, although the 
cutoff energy and spectral index seem to differ somewhat. 
Thus, the Crab Nebula has been well studied by various 
techniques and has been used as a standard calibration 
source among the ground-based 7-ray experiments in the 
TeV region. 

The energy of a primary cosmic-ray particle is es- 
timated by observing the number of secondary parti- 
cles in an air shower experiment as well as the num- 
ber of Cherenkov photons for IACTs, and then by 
comparing these values against the results of detailed 
Monte Carlo simulations, including the detector struc- 
ture and response. The most conventional method for 
estimating the absolute energy scale of primary particles 
may be to compare the flux values between the direct 
(satellite/balloon-borne) and indirect observations. The 
cosmic-ray flux, however, depends inevitably on the de- 
tection technique, analysis method, detector simulation, 
and so on. Therefore, it is very important to develop a 
new method for directly estimating the absolute energy 
scale in a n air shower experiment at TeV energies. 

Clark (1957) anticipated in 1957 that the sun and 
the moon, each with a finite size of 0?5 in diameter, 



cast shadows in the high-energy cosmic-ray flux, respec- 
tively. Actually, the shadowing effect of the sun and 
the moon (hereafter, we call these the sun's shadow 
and the moon's shadow, respectively ) was observed by 
air shower experiments in the 1990s (|Alexand reas 1991; 
lAmenomori et al7l 993). and the sharpness of the ob- 
served shadows was used to estimate their angular res- 
olutions experimentally. In particular, the Tibet air 
shower array, with its high event trigger rate and good 
angular resolution, enables us to use the geomagnetic 
field as a magnetic spectrometer for primary cosmic 
rays at multi-TeV energies. As almost all primary cos- 
mic rays are positively charged, they are bent eastward 
by the geomagnetic field at Yangbajing in Tibet. The 
moon's shadow should then be observed in the west of 
the moon's apparent position, although the position of 
the shadow in relatio n to that of the moon dep ends on 
the cosmic-ray energy ([Amenomori et al. 2000al) . Hence, 
the position and the shape of the moon's shadow al- 
low us to estimate the possible systematic error in 
the absolute energy scale of observed showers. Un- 
til now, the pointing accuracy and angular resolution 
of the Tibet-Ill array have been checked by monitor- 
ing the moon's shadow continuously month by month 
(| Amenomori et al. 20031 ). It is also worthwhile to note 
that the moon's and sun's shadows provide information 
about the cosmic-ray p/p flux ratio (jAchard et al. 20051 ; 
lAmenomori et al. 2007af) and the global structure of the 
interpl anetary magnetic field between the sun and the 
earth (| Amenomori et al. 19931: lAmenomori et al. 1994t 
lAmenomori et al. 1999U lAmenomori et al. 2000aD . 

In this paper, we first discuss the systematic uncer- 
tainties of the Tibet-Ill array and a new method for cal- 
ibrating the absolute energy of primary particles in the 
multi-TeV energy region using the moon's shadow data. 
Based on the results, we report on the 7-ray observation 
of the Crab Nebula and search for pulsed 7-ray emissions 
from the Crab pulsar using the dataset obtained by the 
Tibet-Ill array. 

2. EXPERIMENT 

2.1. Tibet-Ill Air Shower Array 

The Tibet-Ill array shown in Figure Q] was com- 
pleted on the basis of the success of the Tibet- 
I, II and II/HD experi m ents ([Amenomori et al. 19921 ; 
lAmenomori et al. 1999al : lAmenomori et al. 2000bf ) m 
the late fall of 1999. This array consists of 533 scin- 
tillation detectors of 0.5 m 2 , and the detectors on the 
inner side of the array are placed on a lattice with 
7.5 m spacing, covering 22,050 m 2 . A lead plate 0.5 
cm thick is placed on top of each detector to improve 
the angular resolution. Using this array, we succeeded 
in observing 7-ray flares from Mrk 421 and found a 
correlation between Te V 7-ray and X-ray intensities 
(|Amenomori et al. 20031 ). In 2002 and 2003, the in- 
side area of the Tibet-Ill array was further enlarged to 
36,900 m 2 by adding 256 detectors. This full Tibet-Ill 
array has been successfully operating since 2003. In this 
paper, to keep the form of the data the same through- 
out the observation period from 1999 to 2005, we recon- 
structed air shower data obtained from the detector con- 
figuration shown in Figure [T] even for the full Tibet-Ill 
array. 
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Fig. 1. — Schematic view of the Tibet-Ill array operating at 
Yangbajing. Open squares: FT detectors equipped with a fast- 
timing (FT) photomultiplier tube (PMT); filled squares: FT de- 
tectors with a wide dynamic-range PMT; filled circles: density 
detectors with a wide-dynamic range PMT. We have selected air 
shower events whose cores are located within the detector matrix 
enclosed by the dotted line. 

2.2. Event Reconstruction 

The raw data obtained from the Tibet-Ill array system 
mainly consist of the following: a trigger time stamp for 
each event from a global positioning system (GPS) clock 
supplemented by a computer clock; timing and charge 
information from each hit PMT, digitized by a time-to- 
digital converter (TDC) and a charge-sensitive analog- 
to-digital converter (ADC); and calibration data taken 
every 20 minutes. The ADC and TDC counts are then 
converted to the number of particles and the relative tim- 
ing for each detector, respectively, using the calibration 
data. An air shower event is reconstructed as follows. 

The core position of each air shower is estimated using 
the lateral distribution of the number of shower particles 
observed in the array. The density-weighted position of 
the air shower core on the surface of the Tibet-Ill array 

is calculated as (A corc , Fcorc) = f ^r^r > ^fj^rT 
where x$ and j/j are the coordinates of the i-th detector 
and pi is the number density (m~ 2 ) of detected particles. 
In this analysis, we regard the shower size ^ pft as the 
primary energy reference, where size ^2 Pft is defined as 
the sum of the number of particles per m 2 for each FT 
detector. For 7-ray-induced air showers, overall core po- 
sition resolutions are then estimated as 12 m and 4 m at 
median for Pft < 100 and pft > 100, respectively. 

The arrival direction of each shower is estimated as- 
suming that the front of the air shower is conical shape. 
The apex of a cone is then taken to be the estimated core 
position (A corc , y corc ). The average delay time T (ns) of 
shower particles is expressed as a function of the distance 
R (m) from the core position as T = 0.075i?, which is 
optimized by simulations. This gives a cone slope of 1?3 
with respect to the plane perpendicular to the air shower 
arrival direction. 



2.3. Event Selection 

An event trigger signal is issued when an any-fourfold 
coincidence appears in the FT detectors that have each 
recorded more than 0.6 particles within a coincidence 
gate width of 600 ns, resulting in a trigger rate of about 
680 Hz. We collected 2.0xl0 10 events during 1318.9 live 
days from 1999 November 18 through 2005 November 
15 after some quality cuts and event selection based on 
three simple criteria: (1) each shower event should fire 
four or more FT detectors that have each recorded 1.25 
or more particles; (2) among the 9 hottest FT detec- 
tors in each event, 8 should be contained in the fiducial 
area enclosed by the dotted line in Figure [1] If fewer 
than 9 detectors have been hit, they should all be con- 
tained in the fiducial area; and (3) the zenith angle of 
the event arrival direction should be less than 40° . After 
these criteria have been met, the overall angular reso- 
lution and the modal energy of air shower events, thus 
obtained, ar e better than 1 degree and about 3 TeV, 
respectively (|Amenomori et al. 2 003). thereby covering 
the upper portion of the energies measured by IACTs. 

3. MOON'S SHADOW AND PERFORMANCE OF THE 
TIBET-III AIR SHOWER ARRAY 

3.1. Analysis 

Using the dataset described in £|2.3[ we further 
select the events within a circle of the radius 5° 
centered at the moon; this circle is defined as 
the on-source field. An equatorial coordinate sys- 
tem is defined, fixing the origin of coordinates at 
the moon's center. To estimate the background 
against deficits in the moo n's shadow, we adopt 
the equi-zenith angle m ethod (|Amenomori et al. 20031: 
lAmenomori et al. 2 005). which is also used for the 
Crab Nebula observation described in Eight off- 

source fields are symmetrically aligned on both sides of 
the on-source field, at the same zenith angle. In order 
to avoid deficit events that are affected by background 
event contamination of the on-source field, the nearest 
two off-source fields are each set at an angular distance 
9?6 from the on-source field. Other off-source fields are 
located every 3? 2 from the nearest off-source fields. The 
position of each observed event in an on- / off-source field 
is then specified by the angular distance 8 and the posi- 
tion angle <$>, where 6 and <fr are measured from the cen- 
ter and from the north direction, respectively. Using 9 
and (f>, the on-/off-source fields are meshed by 0?05x0?05 
cells, and we count the number of events in each cell. To 
maximize the S/N ratio, we group the cells into new 
on-/off-source bins according to the angular resolution, 
which depends on the shower size pft- The angular 
resolution becomes worst at a threshold energy of around 
1 TeV. This value is, however, smaller than half of the 
angular distance between two adjoining off-source bins, 
i.e., ~ 1?6, so that off-source bins never overlap mutually 
in this background analysis. 

We calculate the statistical significance o f deficits 
or signals using the formula (|Li fc Ma 1983( 1 (N ON - 

eVoFF)/\/ e (^ON + V ff), where N ON and AToff are 
the number of events in the on-source bin and the num- 
ber of background events summed over 8 off-source bins, 
respectively, and e is the ratio of the on-source solid an- 
gle area to the off-source solid angle area (e = 1/8 in this 
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East- Angular Distance -West 

Fig. 2. — Significance map of the deficit event densities observed 
by the Tibet-Ill array for 1318.9 live days, made using the events 
with X)PFT > 10 125 (>~2 TeV), in the square area of 6°x6° 
whose origin is at the apparent center of the moon. The scale at 
right shows the level of significance of the deficit event density in 
terms of the standard deviation a. 

work) . 

In order to investigate the energy dependence, 
the shower size ^ pft is divided by 1 /4 decades 
in 10<^pft <100, and by 1/3 decades in 
100<^ pft<1000, where the lowest air shower size bin 
is omitted from the analysis, because it is close to the 
energy threshold of the Tibet-Ill array and the trigger 
efficiency is estimated to be very low (<1%). Hereafter, 
this partition is commonly used in observations of both 
the moon's shadow and the Crab Nebula. 

Figure [2] shows the experimental significance map of 
the deficit event densities observed with the Tibet-Ill 
array for 1318.9 live days. This map is smoothed using 
the events with Pft>10 125 (>^2 TeV) within a circle 
of radius 0?9, corresponding to the overall angular reso- 
lution for these events. The maximum deficit reaches the 
significance level of 45 a at the center. It is seen that the 
center of the observed moon's shadow is shifted to the 
west by about 0?2 due to the effect of the geomagnetic 
field. 

3.2. Monte Carlo Simulation of the Moon's Shadow 

We have performed a detailed Monte Carlo (MC) 
simulation of the moon's shadow. For the geo- 
magnetic field, we adopt the International Geomag- 
netic Reference Field (IGRF) 9th generation model 
(Macmilla n et al. 20 03) at an altitude < 600 km, and 
connect it to the dipole moment model at an altitude > 
600 km. For the primary particles, we use the chemi- 
cal compositi on obtained mainly from the data of direct 
observations jA sakimo ri et al. 19 98: Sanuki et al. 20001 
Apanascn ko et al. 20011 : lAmenomori et al. 2008T ) in the 
energy range of 0.3 TeV to 1000 TeV. The minimum 
energy of primary particles is set to 0.3 TeV, which is 
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Fig. 3. — Deficit event density map obtained by the MC simu- 
lation. The events with J2 Pft > 10 125 (>~2 TeV) are plotted in 
the square area of 6° x6°. The map is made in the same way as in 
Figure [2] and the scale at right represents the deficit event density 
(degree -2 ). 

low enough to cover the threshold energy of our trigger 
condition. Air shower events are generated at the top 
of the atmosphere along the m oon's orbit aroun d the 
earth, using the CORSIKA code (|Heck et al. 1998f) with 
QGSJET or SIBYLL interaction models. The air shower 
core of each simulated event is uniformly distributed over 
a circular region with a 300 m radius centered at the ar- 
ray; this circle sufficiently covers the area where cosmic- 
ray events are actually triggered in our array. Air shower 
particles generated by primary particles in the atmo- 
sphere are traced until their energies reach 1 MeV. In 
order to treat the MC events in the same way as the 
events in the experimental data, these simulated events 
are distributed among the detectors in the same detec- 
tor co nfiguration as in the Tibet-Ill arra y by the Epics 
code (jAmenomori et al. 20081 : [Kasahara), and are con- 
verted to the same format as the experimental dataset, 
such as the ADC and TDC values at each detector. Af- 
ter air shower reconstruction analysis and data selection, 
we assign the opposite charge to the remaining primary 
particles. These anti-particles are randomly shot back 
toward directions within 20°x20° centered at the moon 
from the first interaction point of the air shower. Here- 
after, we call this the initial shooting direction. The par- 
ticle track influenced by the geomagnetic field between 
the earth and the moon is calculated by the Runge-Kutta 
method of order 4 based on the Lorentz force. If the pri- 
mary particle hits the moon, its initial shooting direction 
should be equivalent to the observed particle direction 
shielded by the moon. Otherwise, it is shot back in an- 
other direction and the particle track is calculated again. 
This routine continues until almost all particles have hit 
the moon. Finally, these initial shooting directions are 
smeared by the angular resolution event by event. In 
this way, the expected moon's shadow is equivalent to 
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the observed moon's shadow. Figure [3] shows the event 
map of the moon's shadow calculated by the MC simula- 
tion for events with X>ft>10 125 (>~3 TeV). This map 
is smoothed using the events within a circle of radius 0?9, 
corresponding to the overall angular resolution for events 
with /?ft>10 1 ' 25 . The MC simulation well reproduces 
the observed moon's shadow, as shown in Figure [21 

3.3. Shape of the Moon's Shadow and Performance of 
the Tibet-Ill Array 

The filled circles in Figure [4] (a)-(f) show the observed 
deficit counts around the moon projected onto the east- 
west axis for each ^ p-pT bin, where Y Pft is the shower 
size defined as the sum of the number of particles per 
m 2 for each FT detector. The representative cosmic-ray 
energy in each ^ p-pT bin is defined as be the logarithmic 
mean of the energy E divided by the atomic number Z, 
assuming the cosmic-ray composition spectrum described 
in §3.21 In these figures, it is seen that the peak posi- 
tion of the deficit counts gradually shifts to the west as 
primary energy decreases due to the influence of the geo- 
magnetic field. Also, the deficit counts become narrower 
as primary energy increases, since the angular resolution 
increases roughly in inverse proportion to Pft- The 
MC results (histograms) are compared with the experi- 
mental data at various size intervals in Figure [4J and are 
in good agreement with the experimental data. 

In order to estimate the peak position of the observed 
moon's shadow, we use the shadow shape obtained by 
the MC simulation. We first express the MC shadows 
shown in Figure|4]with the superposition of two Gaussian 
functions using the least \ 2 method as, 

/mc(0) = Gi(6»;ai,TOi,<7i) 

+G ! 2(6';a2,m2,(T2), (1) 

where Gi(9; dj, m^, cr^) = aje _ ( e_m *) l ai and is the an- 
gular distance from the moon. Here, oi, 02, mi, m.2, 
(Ti and (T2 are the fitting parameters denoting the ampli- 
tudes, means and one-standard deviations of the double 
Gaussian, respectively. It is found that the shadows (e) 
and (f) can be expressed by a single Gaussian, while the 
others are expressed by a double Gaussian. 

Using these fitting parameters, we then estimate the 
peak position of the observed shadow as follows. Keeping 
the form of function /mc(#), we express the observed 
shadow by the equation, 

/Data(#) = G\{9;A\,M\,a\) 
+G 2 (6>; Ai x (02/01), Mi + (m 2 - mi), <7 2 ), (2) 

where A\ and M\ denote the amplitude and mean of the 
Gaussian, respectively, and are free parameters, while 
the others are the coefficients calculated by fitting equa- 
tion (fT]). Using equation ([2]), we can estimate the peak 
position of the moon's shadow. 

We also confirm that the east-west component of the 
geomagnetic field strength is negligible in the part of 
the sky where the moon is visible by the Tibet-Ill ar- 
ray. This means that the north-south displacement of 
the moon's shadow observed by the Tibet-Ill array does 
not depend on the geomagnetic field. The displacement 
of the peak of the moon's shadow in the north-south di- 
rection then enables us to estimate the magnitude of the 
array's systematic pointing error. Figure [5] shows the 



energy dependence of the displacement of the moon's 
shadow in the north-south direction. The filled circles 
denote the experimental data, and the open squares are 
the MC results. The MC simulation well reproduces 
the experimental data. A \ 2 fitting to the data gives 
0?008 ± 0?011 assuming a constant function independent 
of energy. From this, the systematic pointing error is 
estimated to be smaller than 0?011. 

3.4. Calibration of Primary Particle Energy and 
Systematic Errors 

Figure [6] shows the shower size dependence of the dis- 
placement of the moon's shadow in the east-west direc- 
tion, obtained by fitting Figure [4] In this figure, the 
open squares show the MC results using the QGSJET 
model, and are quite consistent with the experimental 
data. The upper scale indicates the logarithmic mean of 
the energy E divided by the atomic number Z (TeV/Z), 
i.e., < \og(E/Z) >, in each Pft bin. One can see 
that the position of the moon's shadow gradually shifts 
to the west as the primary energy decreases due to the 
influence of the geomagnetic field. Hence, the absolute 
energy scale of cosmic rays observed by the Tibet-Ill ar- 
ray can be directly checked by using the geomagnetic 
field as a magnetic spectrometer, as we now discuss. 

First, the MC simulation points are fitted by the func- 
tion k(Y, Pft/100) a to define a standard curvature func- 
tion, resulting in k = —0.183 and A = —0.720, as shown 
by a solid curve in Figure [H where the MC statistical er- 
rors are negligible compared with the experimental data. 

Second, the experimental data (filled circles) are fitted 
by this standard curvature function with the Y Pft shift 
term 

- 0.183 [(1 - AR S )(J2 Pft/100)]- 720 , (3) 

to estimate the possible shift in the Y Pft between the 
experimental data and the MC simulation, as shown by 
the solid curve in Figure[5J where AR$ is the ^ Pft shift 
ratio, resulting in Ai? s = (-4.9 ±9.5)%. We should then 
convert Ai?s to the energy shift ratio ARp as a final re- 
sult. To determine the relationship between ARg and 
AJ?e, and to confirm that this method is sensitive to 
energies, we prepare six kinds of MC event samples in 
which the energy of the primary particles is systemati- 
cally shifted event by event in the moon's shadow sim- 
ulation. These six Ai? E s are ±20%, ±15% and ±8%, 
respectively. In each MC event sample, the Y Pft de- 
pendence of the displacement of the moon's shadow is 
calculated in the same way, and the Y, Pft shift ratio 
Ai?s is estimated by fitting the data to equation |3j) . Fi- 
nally, we get the relation AR E = (-0.91 ± 0.05) AR S 
assuming a linear function. Hence, the systematic error 
in the absolute energy scale ARp with statistical error 
er s tat is estimated to be (±4.5 ± 8.6 sta t)%- 

Furthermore, we investigate two kinds of systematic 
uncertainties with the proposed method. One is that 
the position of the moon's shadow by the MC simulation 
depends on the assumed primary cosmic-ray composi- 
tion. In this simulation, the chemical composition ra- 
tio of primary cosmic rays is estimated based mainly on 
the data obtained by direct observations. These datasets 
should also have statistical and systematic errors. The 
position of the moon's shadow is dominated by the light 
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Fig. 4. — Filled circles show experimental data for deficit counts around the moon projected to the east-west axis for each 
E Pft bin. We use the events contained in the angular band, centered at and parallel to the east- west axis, comparable to the 
E PFT-dependent angular resolution: (a): ±1?4 for 2.94 TeV/Z in 10 1 ' 25 < E Pft < 10 150 ; (b): ±1?0 for 4.20 TeV/Z in 
10 1.50 < J2p FT < 10 175 ; (c): ±0?7 for 6.46 TeV/Z in 10 175 < Epft < 10 2 00 ; (d): ±0?5 for 11.4 TeV/Z in 10 2 00 < Epft < 
10 2 33 ; (e): ±0?3 for 21.6 TeV/Z in 10 2 33 < Epft < 10 267 ; (f): ±0?2 for 45.4 TeV/Z in 10 267 < E Pft < 10 3 00 . The 
solid histograms deno t e the moon's shadow simulation assuming the primary cosmic-ra y composition based on direct observations 
HAsakimori et al. 19981 ; ISanuki et al. 20001 ; lApanasenko et al. 20011; I Amenomori et al. 20081) . 
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Fig. 5. — Dependence of shower size on the displacement of the 
moon's shadow in the north-south direction. The filled circles and 
open squares represent experimental data and the MC simulation, 
respectively. The solid line denotes the fitting to the experimental 
data assuming a constant function, resulting in 0?008±0?011. The 
upper scale indicates the logarithmic mean of E/Z (TeV/Z) in each 
Epft bin. 



Fig. 6. — Dependence of shower size on the displacement of 
the moon's shadow in the east-west direction. The filled circles 
show the experimental data, and open squares represent the MC 
simulation. The solid curve is fitted to the MC events, and dashed 
curves show a ±10% deviation from the solid curve, respectively. 
The upper scale indicates the logarithmic mean of E/Z (TeV / Z) 
in each E PFT bin. 
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Fig. 7. — Dependence of shower size on the displacement of the 
moon's shadow in the east-west direction by the MC simulation 
for the different primary composition models. The filled circles 
and open squares show the experimental data and MC simulation, 
respectively. The solid curve denotes the best-fit curve for the same 
standard composition ratio as in Fig. \E\ (65% P after triggering by 
the Tibet-Ill array, where P means protons). The dashed curves 
show a 6% shift, corresponding to <T sys ti, from the solid curve (see 
text). The downward and upward triangles are the simulated re- 
sults for the proton-rich model (P:75%), and for the heavy-rich 
model (P:55%), respectively. 

component, so that the proton ratio is artificially varied 
by ±10% from a standard chemical composition without 
changing their spectral index, while the other compo- 
nents heavier than helium are varied by +10% in total. 
Figure [7J shows the results for the composition depen- 
dence of primary cosmic rays. The downward triangles 
are the results obtained by the proton-rich model (75% 
protons after triggering by the Tibet-Ill array), while the 
upward triangles are the ones for the heavy-rich model 
(P:55%). These models are fitted by equation ([3]). We 
then obtain er sys ti = ±6% for the systematic error due 
to the difference in chemical composition, as shown by 
the dashed curves in Figure [7J Another systematic un- 
certainty is caused by the difference between hadronic 
interaction models. Figure [8] compares the results for the 
hadronic interaction model dependence by QGSJET with 
those obtained by SIBYLL. It is found that the results by 
the SIBYLL model can be well fitted by equation © ob- 
tained using the QGSJET model. We then obtain <7 sys t2 
= 6% difference between two models. Finally, the dif- 
ference in the energy dependence of the moon's shadow 
between the experimental data and the MC events is 
estimated to be +4.5% (±8.6 sta t ± 6 syst i ± 6/2 syst 2)%. 
This value is within the statistical and systematic er- 
rors. Hence, the absolute energy scale error in the 
Tibet-Ill array is estimated to be smaller than 12% 

= \J&Rl + oftat + °"s ys ti + (o-sy S t2/2) 2 in total averaged 
from 3 to 45 (TcV/Z). 

3.5. On the Energy Estimation of j-Ray Showers 

We established a new calibration method of the abso- 
lute energy scale for cosmic rays based on the moon's 
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Fig. 8. — Dependence of shower size on the displacement of the 
moon's shadow in the east-west direction by the MC simulation 
for different hadronic interaction models. The open squares and 
cross marks are the results obtained by the QGSJET and SIBYLL 
models, respectively. The solid and dashed curves are the best-fit 
results assuming the QGSJET and SIBYLL models, respectively. 

shadow analysis as described above. The air shower in- 
duced by the primary cosmic ray consists of high-energy 
hadronic and electromagnetic cascades. Although several 
plausible hadronic interaction models are prepared in the 
MC simulation, there still remain dependence between 
these models. Therefore, the energy reconstruction from 
the air-shower size depends on hadronic interaction mod- 
els and also the primary chemical composition models. 
These systematic errors were adequately taken into ac- 
count in the absolute energy scale error in this paper. On 
the other hand, the air shower induced by the primary 
7-ray is predominated by the theoretically well-known 
electromagnetic cascades, because the photon cross sec- 
tion for hadronic interactions is approximately two or- 
ders of magnitude smaller than that for the pair creation 
process. Hence, we naturally expect that the absolute 
energy scale error for 7-rays is smaller than 12% which 
is deduced from the moon's shadow observed in the cos- 
mic rays described in i j3.4l In the next section, we will 
provide reliable results on the multi-TeV 7-ray observa- 
tion from the Crab Nebula with the Tibet-Ill air shower 
array finely tuned by the cosmic-ray moon's shadow. 

4. MULTI-TEV 7-RAY OBSERVATION FROM THE CRAB 
4.1. Analysis 

In order to extract an excess of TeV 7-ray events 
coming from the direction of the Crab Nebula (a = 
83?63, S — 22?02), we adopt the same method as 
used in the Mrk 421 analysis in our previous work 
(jAmenomori et al. 2003f ). We call it the equi-zenith an- 
gle method, which is used also for the moon's shadow 
analysis described in i j3.ll The background is estimated 
by the number of events averaged over eight off-source 
bins with the same angular radius as the on-source bin, 
at the same zenith angle, recorded at the same time in- 
tervals as the on-source bin. The nearest two off-source 
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9. — Number of observed air shower events with PFT > 
lO 1,255 (> ~1 TeV) after event reduction for the observation time 
of 1318.9 detector live days as a function of angular distance from 
the Crab Nebula in the azimuthal direction. 
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Fig. 10. — Contour map of significance distribution around the 
Crab Nebula (a = 83?63, 8 = 22?02) for events with X>FT > 
10 1,25 (> ~1 TeV). A clear peak excess is seen at the center po- 
sition Aacos(A<5) = A<5 = 0°, where Aa and A<5 are the relative 
right accension and declination, respectively, from the Crab Neb- 
ula. The cross mark indicates the pointing error by a point spread 
function fitting. 

bins are set at an angular distance 6? 4 from the on-source 
bin to avoid a possible signal tail leaking into these off- 
source bins. Other off-source bins are located every 3?2 
step from the nearest off-source bins. The search window 
radius is expressed by 6.9/ VzT/°ft degrees as a function 
of ^2 Pft i which is shown to m aximize the S/N ratio b y 
MC study of 7-ray observation (jAmenomori et al~2 003) . 

The number of events after the event reduction is plot- 
ted in Figure [9] as a function of angular distance from the 
Crab Nebula in the azimuthal direction. A clear peak of 
7-ray signals from the Crab Nebula is seen at 6.3 a statis- 
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Fig. 11. — Scatter plot of the shower size Pft an d t ne energy 
of primary 7-rays, where a differential power-law spectrum of the 
form E~ 2,6 starting at 0.3 TeV is assumed for primary 7-rays. For 
details, see text. 

tical significance above the flat cosmic-ray background. 
Figure [10] is a significance map around the Crab Neb- 
ula. The peak excess is seen at the Crab Nebula posi- 
tion. The pointing error as shown by a cross mark in 
Figure [TU] is estimated to be Aa = +0.13 ± 0.08 and 
AS = +0.01 ± 0.09 by the point spread function fitting. 
This is consistent with the Crab's position within sta- 
tistical error. As difference between 7-ray induced air 
showers and cosmic-ray induced ones does not affect the 
pointing accuracy essentially, we estimate our pointing 
accuracy to to be 0?011 deduced from the moon's shadow 
analysis described in £13.31 

4.2. Monte Carlo Simulation of "f-Ray Observation 
from the Crab Nebula 

The performance of the Tibet-Ill array has been stud- 
ied by a full MC simulation using the CORSIKA code 
(|Heck et al. 19981 ) for even t generati on in the atmo- 
sphere and the Epics code llKasaharaD for the respons e 
of the scintillation detector (|Amenomori et al. 20031) . 
These procedures are essentially the same as in the case 
for the moon's shadow described in ^3.21 In the simu- 
lation for 7-ray observation from the Crab Nebula, pri- 
mary 7-rays, assuming the energy spectrum of a power- 
law type in the energy region of 0.3 TeV to 1000 TeV, are 
thrown along the diurnal motion of the Crab Nebula in 
the sky. The air shower events generated are uniformly 
distributed over circle with a 300 m radius centered at 
the Tibet-Ill array. Shown in Figure [TT] is the scatter 
plot of shower size ^2 Pft and the energy of 7-rays com- 
ing from the Crab direction. The filled circles and error 
bars stand for the logarithmic mean of the energy and 
one-standard deviation of the logarithmic Gaussian, re- 
spectively. The one-event energy resolution is estimated 
to be approximately (—40/ + 70)% at 10 TeV, and ap- 
proximately ±100% in the region of a few TeV. 

The Crab Nebula can be treated as a point-like source 
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Fig. 12. — Distribution of excesses as a function of the opening 
angle relative to the Crab Nebula direction 9. The filled circles 
and shaded histograms stand for the experimental data and the 
MC events with £ p PT > 10 125 (> ~1 TeV), respectively. 



at the TeV energy region. To investigate the point spread 
function of the Tibet-Ill array, we compared the 9 dis- 
tribution of the Crab Nebula between the experimental 
data and the MC events, where 8 is the opening angle 
relative to the Crab Nebula direction. Figure [12] shows 
the distribution of the excess events as a function of 9 
for events with pft > 10 1 ' 25 . The experimental data 
agree well with the MC simulation assuming the point- 
like source. 

4.3. Energy Spectrum of 7- Rays from the Crab Nebula 

The 7-ray flux from the Crab Nebula is estimated by 
assuming a power-law spectrum f(E) — aE^ . The best- 
fit values cto and /3q are given by minimizing a x 2 func- 
tion, changing a and (3: 



6 

E 



N° bs - Nf m (a,P) 



_obs 



(4) 



where N° , a° and Nf lm (a,f3) are the observed num- 
ber of excess counts, its error and the number of remain- 
ing MC events after the analysis assuming the spectrum 
f(E) = ctE$ , respectively, in the z-th X] Pft bin among 
the six pft bins between 10 125 and 10 3 00 defined in 
N3- 11 In order to estimate N!* lm (a, 0) in the same way as 
experimental data, simulated secondary particles are in- 
putted to the detector response simulation. Then, we ob- 
tain the expected Nf lm (a, /3) for the z-th J2 Pft bin after 
the event reconstruction and event selections in the same 
way as experimental data. Here, the expected N!* lm (a, (3) 
includes the energy resolution effect by the detector re- 
sponse simulation. 
Subsequently, the differential 7-ray flux for the i-th 
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Fig. 13. — Differential energy spectrum of 7-rays from 
the Crab Nebula obtained using the data collected from 
1999 November to 2005 November with the Tibet-Ill ar- 
ray in compariso n with the results from IACTs : Whipple 
HHillas et al. 19981). HEGRA IIAharonian et al. 2004). CANGA- 
ROO III ( Enomo toet al. 2009) H.E.S.S. llAharonian et al. 20061) 
and MAGIC l|Albert et al. 20081 ). The Tibet-Ill upper limit is 
given at t he 90% conf idence level, according to a statistical pre- 
scription (Hclcnc 1983). 

Y^, Pft bin is calculated by the following equation: 



fi(Ei) 



n: 



obs 



Vf lm (a ,/?o) 



A a T(ao,A>) 



E 00 dE 



E 



00 



(5) 



where V^jJ^ao, /?o) denotes the total number of MC 
events generated at the top of the atmosphere along one 
diurnal motion assuming the spectrum f(E) = aoE^°, 
J^i m E 13 " dE is the normalization factor of N^lcto, Po), 

-E^"n denotes the minimum energy of simulated air 
shower events (0.3 TeV), S'sim denotes the area of core 
location distribution by the simulation (300 m x 300 m 
x 7r), T b s denotes the observation live time, and Ei de- 
notes the representative energy defined as the logarith- 
mic mean of the energy calculated by the MC simulation 
for the i-th ^ pft bin. 

Figure [TBI shows the differential energy spectrum of the 
Crab Nebula observed by the Tibet-Ill array together 
with the spectra o btained by IA CTs, including Whippl e 
ijHillas et al. 19981), HEGRA ( Aharoni an et al. 20041) , 
CANGARO O HI (lEnofnoto et al. 20061) , 
H.E.S.S. (lAharonian et al. 2006h and MAGIC 
(| Albert et al. 20081 ). The differential flux for each 
^2 Pft bin is presented in Table [1] Finally, this 
energy spectrum is fitted by the least \ 2 method 
assuming f(E) — a(E/3 TeV)' 3 , and then we ob- 
tain the differential power-law spectra as (dJ /dE) = 
(2.09±0.32) x 10- 12 (£V3 TcV)- 2 - 96±014 cm- 2 s- 1 TcV- 1 
in the energy range of 1.7 TeV to 40 TeV. Note that 
the absolute energy scale error in the Tibet-Ill array is 
experimentally estimated to be smaller than ±12% by 
the moon's shadow observation described in ^3.41 The 
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TABLE 1 

Logarithmic mean of energy and differential flux for each ^Zpft bin as shown in Fig. 1131 



Dpft 


Energy 


Non 


cTVoff 


Differential Flux 






[TeV] 






[cm" 2 s" 1 TeV" 1 


] 




1.71 


1935499 


1931547 


(8.72 ±3.25) x 10- 


12 


10 1.50 _ 10 1.75 


2.89 


1382356 


1377139 


(2.70 ±0.643) x 10" 


-12 


10 1.75 _ 10 2.00 


4.84 


444504 


442074 


(5.52 ± 1.60) x 10- 


13 


Xo 2 00 - 10 2 - 33 


8.29 


134509 


133362 


(1.03 ±0.348) x 10" 


-13 


10 2.33 _ 1Q 2.67 


18.6 


21530 


21138 


(1.06 ±0.417) x 10" 


-14 


10 2.67 _ 1Q 3.00 


39.5 


3923 


3844 


(6.64 ±5.57) x 10- 


16 


> 10 3 00 


107 


1558 


1569 


< 1.10 x 10~ 16 





energy scale uncertainty corresponds to (—28/ + 46)% 
in the absolute 7-ray flux, assuming the spectral index 
—2.96, which is our best-fit value. Our energy spectrum 
in this work is consistent with other observations made 
by IACTs, such as HEGRA and H.E.S.S., in the same 
energy range between 1.7 TeV and 40 TeV. 

The previous flux measurement 

(| Amenomori et al. 1999a| ). with the Tibet-HD ar- 
ray of 5,175 m 2 and an effective running time of 502.1 
live days, is approximately double this measurement. 
In order to properly estimate the difference between 
the previous work and the present one, we give a 
re-fit to both data points from 2.8 TeV to 20 TeV in 
the overlapping energy region assuming a power-law 
spectrum. The previous (Tibet-HD) and present 
(Tibet-Ill) energy spectra are expressed as (dJ/dE) = 
(5.04±0.94) x 10 Ll2 (£/3 TcV)- 2 - 85±0 - 20 c m - 2 s- 1 TeV- 1 



and (dJ/dE) 
10- 12 (£/3 TeV)- 3 00±a25 



(2.35 ± 0.49) x 
cm "s -1 TeV -1 , respectively. 
The flux and spectral index differences between them 
are estimated to be (2.69±1.06) x lO^cm^s^TeV -1 
and 0.15±0.32, respectively. As a result, the combined 
st atistical deviation between them is calculated to be 
^(2.69/1. 06) 2 + (0.15/0.32) 2 cr = 2.6cr. Although we 
have updated the MC simulation in this analysis, we 
cannot find any systematics to explain this difference. 
Hence, we conclude that the higher flux observed in 
our previous measurement may have been caused by a 
statistical signal fluctuation. 

4.4. Time Variability 

We divided our dataset from 1999 November to 2005 
November into six phases, as summarized in Tabled to 
examine the time variability of the flux intensity. Each 
phase corresponds to approximately one calendar year. 
We used slightly different calibration parameters for each 
phase, because we usually calibrate the scintillation de- 
tectors of the Tibet air shower array late in the fall of 
every year. Unfortunately, some of the blank periods 
seen in Table [5] mostly coincide with the detector cali- 

TABLE 2 

Definition of six phases from 1999 November to 2005 
November as shown in Fig. O 



Phase 


Period 




Live Time 










[days] 


1 


Nov. 18, 1999 - 


Jun. 


29, 2000 


173.1 


2 


Oct. 28, 2000 - 


Oct. 


11, 2001 


283.7 


3 


Dec. 05, 2001 - 


Sep. 15, 2002 


201.8 


1 


Nov. 18, 2002 - 


Nov. 


18, 2003 


259.1 


5 


Dec. 14, 2003 - 


Oct. 


10, 2004 


123.6 


6 


Oct. 19, 2004 - 


Nov. 


15, 2005 


277.6 



bration periods, periods in which the air shower array 
was upgraded or when the data acquisition system was 
experiencing problems. The upper panel (a) in Figure fT4l 
shows the time variability of the 7-ray fluxes from the 
Crab Nebula at 3 TeV. We found no evidence for the 
time variability of flux intensity from the Crab Nebula, 
as we can give a good \ 2 fit to these fluxes by a con- 
stant function (% 2 /d.o.f. = 6.55/5), where d.o.f. means 
degrees of freedom. In order to check the possible sys- 
tematics, the time variability of the deficit event rates of 
the moon's shadow is also demonstrated as shown by the 
middle panel (b) in Figure Q3] The deficit event rates of 
the moon's shadow from 1999 November to 2005 Novem- 
ber are very stable within a statistical error ±6% year 
by year. A fitting to the daily deficit event rate averaged 
over a phase assuming a constant function is consistent 
with a flat hypothesis (x 2 /d.o.f. = 4.82/5). The lower 
panel (c) in Figure [14] shows the time variability of the 
north-south displacement of the moon's shadow, which 
is a reference to the absolute pointing error described 
in §3.21 It is also very stable within ±0?04 during our 
observation period, and is consistent with a flat hypoth- 
esis (x 1 /d.o.f . — 2.09/5). These systematics, estimated 
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Fig. 14. — Time variability of the Crab Nebula and the moon's 
shadow observed by the Tibet-Ill array with Pft > 10 1 ' 25 be- 
tween 1999 November and 2005 November, (a): Differential flux of 
the Crab Nebula at 3 TeV. (b): Daily deficit event rate averaged 
over one phase of the moon's shadow, (c): North-south displace- 
ment of the moon's shadow. 
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from observations of the moon's shadow, are obviously 
negligible in the Crab Nebula observation. 

4.5. Search for j-Rays from the Crab Pulsar 

The rotation period of the Crab pulsar is 33 ms, as 
inferred from radio, optical and X-ray observations. A 
pulsed emission with that rotation period at the GeV 
energy region has b een detected by EGRET on board 
the CGRO satellite ()Fierro et al. 19981 ). whereas several 
observations have reported no e vidence for pulsed 
emissions greater than 10 Ge V dLessard et al. 2000 ; 

De Naurois et al. 20021 lAharonian et al. 2004 ; 

Aharonian et al. 20071: I Albert et al. 2008fh The 
emission models of high-en ergy pulsed 7-rays a re mostly 
based on t he outer gap (ICheng et al. 19861) and the 
polar cap (jDaughertv fc Harding 19821 ) models. The 
upscattered pulsed 7-ray flux is also calculated by 
the inverse-Compton process and the photon-photon 
absorption process assuming infrared photon field 
models. A model predicts an excessive flux around 1 ~ 
10 TeV, depending on the infrared photon field models 
(jHirotani fc Shibata 200ll ). Here, we present a search 
for pulsed 7-rays from the Crab pulsar at energies from 
a few TeV to 100 TeV using the Tibet-Ill array. 

The arrival time of each event is recorded using a 
quartz clock synchronized with GPS, which has a pre- 
cision of 1 fj,s. For the timing analysis, all arrival times 
are converted to the solar sy stem barycenter fram e using 
the JPL DE200 ephemeris (jStandish et al. 1982f) . The 
Crab pulsar ephemeris is calculated using the Jodrell 
Bank Crab Pulsar Monthly Ephemeris (jLvne et al. 19931 
iLvne et "ail ). The corrected arrival time of each event 
is calculated to the rotation phase of the Crab pulsar, 
which takes into account of the period derivative P of 
the period P month by month. 

Figure [15] shows the distribution of events for each 
phase in two rotation periods of the Crab pulsar. 
The distribution is consistent with a flat distribution 
(x 2 /d.o.f. = 18.1/19). No significantly pulsed signal is 
found in observations for events with Pft > 10 125 (> 
~1 TeV). The phase analysis is performed for each Pft 
bin to examine the energy depende nce. TableBlshows th e 
statisti cal results by th e Z\ test (jBuccheri et al. 1 983). 
H test ()De Jager 19941) and least \ 2 test- Almost all the 
statistical tests show that the phase distributions are uni- 
form within a 3 a significance level. We estimate the 3 a 
flux upper limit on the pulsed emission from the Crab 
pulsar using the H test (|De Jager 19941) as 

x 3a = (1.5 + 10.75) (0.174if) 017+0 - 14 * 

x exp{(0.08 + 0.15(5) 

x(log 10 (0.174#)) 2 }, (6) 

where S is the duty cycle of the pulsed component, as- 
suming S is 21% for the Crab pulsar. Exposure from the 
Crab pulsar to the Tibet-Ill array is estimated using MC 
simulation, assuming the differential energy spectrum for 
7-ray emission has a spectral index of —2.6. Upper lim- 
its are compared to previous results inferred from other 
experiments, as shown in Figure [TBI 

5. SUMMARY AND PROSPECTS 

We have been successfully operating the Tibet-Ill 
air shower array at Yangbajing in Tibet, China since 
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limits on the pulsed 7-ray flux 
observed by the Tibet-Ill ar- 
solid line) , toge ther with re- 
I ILessard et al. 20001), CELESTE 
lIDe Naurois et al. 2001). HEGRA l l Aharonian et al. 2004TI . 
H.E.S.S. [lAharonian et al. 20071) and MAGIC 

l|Albert et al. 20081V The long-dashed curve and dashed 
curve represent the fluxes expected from the outer gap and polar 
cap models, respectively. 



TABLE 3 



2|- 



AND //-TEST (PROBABILITIES) ARE CALCULATED FOR A FLAT PHASE 
DISTRIBUTION. 



SPFT 



lr JTT25 _ 10 l 
10 1.50 _ 10 1.75 
10 1.75 _ 1Q 2.00 
IQ2.00 _ 10 2 - 33 
10 2.33 _ 10 2.67 
10 2.67 _ 1Q 3 

> 10 3 



Tran 



\ 2 /d.o.f. 



-i3.no 



> 10 



3.00 

T723 - 



0.97 (0.49) 
1.21 (0.24) 
0.81(0.70) 
0.35 (0.96) 
1.41 (0.11) 
0.80 (0.71) 
0.60 (0.91) 



-TrZ- 

£2_ 



H 



9.62 (0.047) 
7.64 (0.11) 
2.54 (0.64) 
2.30 (0.68) 

9.68 (0.046) 
3.67 (0.45) 
1.11 (0.89) 



9.62 (0.021) 
7.64 (0.047) 

4.49 (0.17) 
6.14 (0.086) 
14.56 (0.0030) 
6.09 (0.088) 

4.48 (0.17) 



0.95 (0.52) 8.41 (0.078) 8.87 (0.029) 



12 



Amenomori et al. 



1999. Using the dataset collected by this array 
from 1999 November through 2005 November, we ob- 
tained the differential energy spectrum of 7-rays from 
the Crab Nebula as (dJ/dE) = (2.09 ± 0.32) x 
10~ 12 (£;/3 TeV)- 2 - 96±014 cm- 2 s- 1 TcV _1 in the energy 
range of 1.7 TeV to 40 TeV. This result is consistent 
with data obtained by IACTs, and is statistically consis- 
tent with our previous result within 2.6 a. No evidence 
is found for time variability of flux intensity from the 
Crab Nebula at multi-TcV energies in comparison with 
the long-term stability of the moon's shadow. We also 
searched, unsuccessfully, for pulsed 7-rays from the Crab 
pulsar at multi-TeV energies. 

In this paper, we have carefully analyzed the moon's 
shadow observed with the Tibet-Ill array to calibrate the 
energy of primary cosmic rays directly. In general, this 
energy is indirectly estimated by measuring shower size 
in air shower experiments. The cosmic-ray beams coming 
from the moon's direction are bent by the geomagnetic 
field, so that the moon's shadow should shift to the west 
depending on the primary energy. We tried to directly 
estimate the primary energy by measuring the displace- 
ment of the moon's shadow. This energy scale estima- 
tion is the first attempt and obtained that the systematic 
error in the absolute energy scale observed by the Tibet- 
Ill array is estimated to be less than ±12% at energies 
around 10 TeV. The array's systematic pointing error 
is also estimated to be smaller than 0?011. The long- 
term stability of the deficit rate of the moon's shadow 
was within a statistical error ±6% year by year, thus 
confirming the stability of the array operation. This cal- 
ibration method is very unique and will be important to 
ground-based TeV 7-ray observations. 



In the near future, we will set up a 10,000 m 2 
water- Cherenkov- type muon detector (MD) array in 
the ground beneath the Ti b et air shower (AS) ar- 
ray (lAmenomori et al. 2007bt lAmenomori et al. 2007d ; 
I Amenomori et al. 2007of K This Tibet MD array will 
significantly improve 7-ray sensitivity of the Tibet air 
shower array above 10 TeV by means of 7/hadron sep- 
aration based on counting the number of muons accom- 
panying each air shower. The energy spectrum of the 
Crab Nebula can be surely measured up to several hun- 
dred TeV, if extended, with a low background level, us- 
ing the Tibet AS+MD array. This new array will en- 
able us to survey not only the known sources but also 
new sources in the northern sky above 10 TeV, and may 
be superior to IACTs for observing diffuse 7 -ray sources 
(|Amenomori et al. 20061 : lAbdo et al. 20071) and diffuse 
7-rays from the ga lactic plane (jAmenomori et al. 20021 : 
lAtkins et al. 20051 ). owing to its wide field of view and 
high rejection power for hadronic showers. 



The collaborative experiment of the Tibet Air Shower 
Arrays has been performed under the auspices of the 
Ministry of Science and Technology of China and the 
Ministry of Foreign Affairs of Japan. This work was sup- 
ported in part by a Grant-in- Aid for Scientific Research 
on Priority Areas from the Ministry of Education, Cul- 
ture, Sports, Science and Technology, by Grants-in-Aid 
for Science Research from the Japan Society for the Pro- 
motion of Science in Japan, and by the Grants from the 
National Natural Science Foundation of China and the 
Chinese Academy of Sciences. 
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